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Monotone finite difference schemes for anisotropic diffusion 
problems via nonnegative directional splittings* 

Cuong Ngo* Weizhang Huang* 


Abstract 

Nonnegative directional splittings of anisotropic diffusion operators in the divergence form 
are investigated. Conditions are established for nonnegative directional splittings to hold in a 
neighborhood of an arbitrary interior point. The result is used to construct monotone finite 
difference schemes for the boundary value problem of anisotropic diffusion operators. It is 
shown that such a monotone scheme can be constructed if the underlying diffusion matrix is 
continuous on the closure of the physical domain and symmetric and uniformly positive definite 
on the domain, the mesh spacing is sufficiently small, and the size of finite difference stencil 
is sufficiently large. An upper bound for the stencil size is obtained, which is determined 
completely by the diffusion matrix. Loosely speaking, the more anisotropic the diffusion matrix 
is, the larger stencil is required. An exception is the situation with a strictly diagonally dominant 
diffusion matrix where a three-by-three stencil is sufficient for the construction of a monotone 
finite difference scheme. Numerical examples are presented to illustrate the theoretical findings. 
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1 Introduction 


We consider the finite difference solution of the boundary value problem 

J—V • (DVu) = /, in Q 
\u = g, on 

where O = (0,1) x (0,1), / and g are given functions, and the diffusion matrix, 

|"a(x, y) b(x,y) 

[b{x,y) c[x,y)\ ’ 


(1) 


(2) 


is assumed to be continuous on Q e llU 30 and symmetric and uniformly positive definite on 0. 
We are interested in the case where D depends on spatial location (heterogeneous diffusion) and has 
unequal eigenvalues at least on some portion of O (anisotropic diffusion). For this case, 0 is often 
called a heterogeneous anisotropic diffusion problem. Anisotropic diffusion arises from various areas 
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of science and engineering, including plasma physics [9j, petroleum reservoir simulation EES, and 
image processing mm- 

When a conventional method such as a finite element, a finite difference, or a finite volume 
method, is applied to this problem, spurious oscillations can occur in the numerical solution. This 
difficulty can be overcome by numerical schemes satisfying a discrete maximum principle (DMP). 
The development and studies of DMP-preserving schemes have received considerable attention in 
the past particularly in the context of finite element and finite volume methods, e.g., see mmum 
Emnum H3UHUI5UI6UI7UI8UI9U21U25U26U271 [28U29U30U321 [33] . On the other hand, relatively 
less work has been done in the context of finite difference (FD) methods. The major effort has 
been made to study monotone (or called positive-type or nonnegative-type) schemes, which belong 
to a special class of DMP-preserving schemes with the coefficient matrix (or the Jacobian matrix 
in the nonlinear case) of the corresponding algebraic system being an M- matrix (e.g., see [301] 
for the definition of M-matrices). For example, Motzkin and Wasow [21] prove that a monotone 
FD scheme exists for any linear second-order elliptic problem when the mesh is sufficiently fine. 
Greenspan and Jain [ 8 ], using nonnegative directional splittings (see the definition below), propose 
a way to construct such schemes for elliptic operators in the nondivergence form 

L[u] = a(x,y)u xx + 2b(x,y)u xy + c(x,y)u yy with b(x, y) 2 < a(x, y)c(x, y). (3) 

The results are extended to elliptic problems in the divergence form 0 by Weickert m ■ Consistent 
and stable monotone FD schemes are shown to be convergent (to the solution) for linear second- 
order elliptic problems by Bramble et al. [2] and (to the viscosity solution) for nonlinear second- 
order degenerate elliptic or parabolic partial differential equations by Barles and Souganidis [T] . 
Oberman [23] studies degenerate elliptic schemes (a special type of monotone scheme) for a general 
class of nonlinear degenerate elliptic problems. 

The objective of this work is to study the construction of monotone FD schemes for problems 
in the form of 0 using nonnegative directional splittings. Weickert’s results are improved in two 
aspects. We first present a condition under which nonnegative directional splittings hold for a 
neighborhood of an arbitrary point. As we will see below, it is necessary to consider nonnegative 
directional splittings in a neighborhood for the construction of monotone FD schemes. We then 
extend the result to the situation where the coefficient b(x, y ) can change sign over the domain. To 
be more specific, we recall that Weickert [3lj considers the directional splitting 

V • (BX7u)(x 0 ,y 0 ) = d x ('y 0 d x u)(x 0 ,y 0 ) + dp^idpu^xo^o) + d y (^ 2 d y u){xQ, y 0 ), (4) 

where (xo,yo) is a given point in Q, dp = (sin(/3), cos(/3)) T • V, 7 * = 7 i(x,y) (i = 0 , 1 , 2 ) are 

functions, and /3 is a constant (angle) (but can vary with (xo,yo))- Weickert [3T], Page 90] shows 

that the condition 

a(x 0 , yo) - b(x 0 , y 0 ) cot(/3) > 0 and c(x 0 , y 0 ) - b(x 0 , y 0 ) tan(/3) > 0 (5) 

is sufficient to guarantee that 7 j (i = 0 , 1 , 2 ) are nonnegative at (xo,yo) (nonnegative directional 

splitting). Unfortunately, such a point-wise result is insufficient for the construction of monotone 
FD schemes for the divergence form Q since an FD discretization has to use the information of the 
coefficients 7 * in a neighborhood of any mesh point. To avoid this difficulty, we develop a condition 
under which nonnegative directional splittings hold in a neighborhood of a given point. Moreover, 
we study how the condition can be extended to the situation where b(x, y) changes sign over the 
domain. An upper bound for the stencil size is obtained, which is completely determined by the 
diffusion matrix. Loosely speaking, the more anisotropic the diffusion matrix is, the larger stencil is 
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required. An exception is the situation with a strictly diagonally dominant diffusion matrix where 
a three-by-three stencil is sufficient for the construction of a monotone finite difference scheme (cf. 
Theorem 3.2). 

An outline of this paper is as follows. Nonnegative directional splittings for 0 are studied in 
l|2j Construction of monotone FD schemes using the results in §[2] is discussed in fj3j Numerical 
examples are presented in §|4j followed by the conclusions in ^5} 


2 Nonnegative directional splitting 

In this section we consider nonnegative directional splittings for problems in the form of |l]) . Recall 
that our goal is to use the splittings to construct monotone FD schemes. Thus we are interested in 
the situations where such a splitting holds for a neighborhood of a given point. We first consider the 
case where b(x, y) does not change sign and then extend the result to the sign changing situation. 
Finally, we establish a condition under which the nonnegative directional splitting holds uniformly 
over the domain for any neighborhood of a fixed size. The result is needed in constructing monotone 
FD schemes on a uniform mesh. 

Lemma 2.1. Consider an arbitrary point (a?o,yo) ^ an d a neighborhood Do °.f ( x 'o ; Vo)■ We 
assume that b(x,y) does not change sign on Do- If there exists a constant /3 satisfying 

sup ^ J < tan /3 < inf ^ , for b > 0 and 6^0 on Do 

no a(x,y) - 0 b(x,y)’ 

' SU P 77^4 < tan 0< inf hi pA for b < 0 and b ^ 0 on Q 0 ( 6 ) 

n ,b{x,y)~ H ~ tio a(x,y)’ ~ 

(3 is any finite value with sin(/3) cos(/3) f 0, for 6 = 0 on Do 
then we have the nonnegative directional splitting 


V • (DVu) (x 0 ,y 0 ) = d x ('y 0 d x u)(xo,yo) + ^(71 dpu)(x 0 ,yo) + d y (j 2 d y u)(x 0 ,y 0 ), (7) 


where the coefficients 70 , 71 , and 72 , all nonnegative on Do, are given by 


7o (x,y) 
< li(x,y) 
l2(x,y) 


a(x, y ) — b(x, y) cot (3 > 0, 
b(x, y) 

— v > 0, 

cos j3 sin (3 

c(x, y) — b(x, y) tan (3 > 0. 


( 8 ) 


Proof. For the situation with b(x,y) e 0 on Do, it is obvious that 0 and 0 hold with any 
(3 satisfying sin(/3) cos(/3) f 0. For other situations, condition ([ 6 ]) implies that 0 < /3 < \ or 
— < j3 < 0. Since D (x,y) is symmetric and uniformly positive definite on Do C D, we have 
a(x,y) > 0 and c(x,y) > 0 for each (x,y) 6 Do- For each (x,y) £ Do, we denote the eigenvalues 
of of D by Ai and A 2 (without loss of generality we assume Ai > A 2 ) and the angle formed by 
the principal eigenvector and the x-axis by if. Then, the diffusion matrix at (x, y) has the eigen- 
decomposition as 


— Ai 


cos(if) 
sin (if) 

Ai cos 2 (if) + A 2 sin 2 (V>) 
(Ai — A 2 ) cos(V’) sin(if) 


[cos(V>) sin(V>)] + A 2 


— sin(^) 
cos (if) 

(Ai - A 2 ) cos(if) sin(V’) 

Ai sin 2 (if) + A 2 cos 2 (-i/’) 


[— sin(-i/’) cos(if)] 


(9) 
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In the following, we assume 6 > 0 on ho- (The analysis for the case with b < 0 is similar.) 
In this case, if G [0, j\. We now show that V • (DVu) can be written into the form of (|7j) with 
nonnegative 70 , 71 , and 72 . By direct calculation, we have 

cos 2 (3 
cos (3 sin (3 

'1 O' 

0 0 


7o 


c>s( 7 i dpu) = V • ( 71 
dxho d x u) = V 

d y <Ci 2 d y u ) = V • ( 72 
Then, 0 is equivalent to 

V • (B Vii) = V • 


cos j3 sin f3 
sin 2 f3 


Vu 


Vu 


0 0 
0 1 


Vu 


7 0 + 71 cos 2 /? 71 cos (3 sin f3 

71 cos (3 sin j3 71 sin 2 f3 + 72 


Viz 


which requires 


a b 


7 o + 71 cos 2 /3 71 cos (3 sin f3 

b c 


71 cos (3 sin (3 71 sin 2 f3 + 72 


This can be expressed into 


1 cos 2 /3 0 


7o’ 


a 

0 cos /3 sin j3 0 


7i 

= 

b 

_0 sin 2 /3 1 


.72. 


c 


Since 0 < (3 < |, we can solve the above equation and get ([ 8 ]). The condition (| 6 j) implies that 
70 > 0 and 72 > 0 on Ho- From ([9]), we have 

b = (Ai — A 2 ) cos i/j sin if. 


Thus, 


7i = 

since if G [0, |] and f3 G (0, §). 


b 

sin /3 cos (3 


(Ai — A 2 ) cos if sin if 
cos /? sin (3 


> 0 , 


□ 


The above lemma shows that for the case where b does not change sign, 0 has a nonnegative 
directional splitting on Qq if ([ 6 ]) is satisfied. In the following lemma, we consider the general case 
where b changes sign on Ho- Hereafter, we use the superscripts “ + ” and “ to indicate the regions 
associated with b(x, y) > 0 and b(x, y) < 0, respectively. For example, we denote 

^0 = {( x > v) G ^0: &(z, y) > 0}, fig = {( x >y) e ^0: K x ,y ) < °}- 


Remark 2.1. The set Hq" or Hq may be empty. When this happens, many operations associated 
with such an empty in the following analysis may not make sense. For this reason, we assume that 
the formulas with an empty set will be ignored. This remark also applies to other empty sets. □ 


Lemma 2.2. Consider an arbitrary point (xQ,yo) in II and a neighborhood Ho °f ( x o,yo)- If 
there exist constants j3\ and fa satisfying 


Hx,y) _ a ^ • f c( x ,y) 

sup —-- < tan pi < ml —- 

n+ a{x,y) n+ b{x,y) 

c ( x ,y) / , a ^ f b ( x >y) 

sup — -- < tan fa < mf —--, 

n - b ( x ^y) n- a{x,y) 


( 10 ) 
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then we have the nonnegative directional splitting 

V • (DVu) (x 0 ,y 0 ) = d x (j 0 d x u)(x 0 ,yo) + %(7i’ d/ 3 1 u)(x 0 ,y 0 ) 

+ dfohT d/32 u )( x o,yo) +dy{j 2 d y u){x 0 ,y 0 ), (11) 

where 70 (x,y), / yf(x,y), 7 f(x,y), and 72 (x,y), all nonnegative on Qq, are given by 


7o (x , 

.2/) = ! 

f a{x,y) - b{x 

,y) cot Pi, 

for b(x, y) 

> 

0 

(12) 

{a(x,y ) - b(x 

,y) cot p 2 , 

for b(x, y) 

< 

0 


■ y) = \ 

[ b(x,y) 

for 

b(x, y) > 0 



(13) 

7i~ (x, 

cos Pi sin Pi' 






lo, 

for 

b(x, y) < 0 






\o, 

for 

b(x,y) > 0 




7i~ (*. 

-y) = \ 

1 b(x,y) 





(14) 


1 

l cos p 2 sin p 2 : 

J U1 

y) ^ u 




12{X. 

■ y) = | 

\c{x,y) - b(x. 

y) tan Pi , 

for b(x, y) 

> 

0 

(15) 


{c{x,y) - b(x. 

y) tan , 

for b(x, y) 

< 

0 . 


Proof. For any given constants e > 0 and M > 0 we define a decompose of the diffusion matrix as 


where 


B(x,y) = Di(x,y) +B 2 {x,y), 



bi 

ci 


O2 


a 2 62 
b 2 c 2 


a — e b + e/M 
b + e/M c — e 

e 0 
0 e ’ 

e -e/M 
-e/M e 

a — e b 


for b(x, y) > 0 

for b(x, y) < 0 
for b(x, y) > 0 

for b(x, y) < 0. 


(16) 


(17) 


Then we have 

Since 


V • (DVu) = V • (DiVrt) + V • (ID^Vit). 

b\ > 0 and b 2 < 0, V(x,y) € Hq 


(18) 


we can apply Lemma 2.1 to each term on the right-hand side of (18). To this end, we notice that 
(10) implies, for sufficiently small e and large M, that 


b b + e/M . p c e c 

sup - < sup- < tan pi < ml -—— < mf -, 


HI 


a 


+ a — e 


n. 


+ b + e/M 




(19) 
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c C-e . - 6 . , & 

sup - < sup —-— < tanP2 < int- < ini 


n, 


- b 


From (16), we have 


b\ h b + e/M 

sup — = 0 < sup — = sup 


Q- 


O"I- ^1 o 


+ a — e 


Q“ a — e n~ a 


inf ^ = +cx) > inf ^ = inf —-—. 
n- h n+h n+b + e/M 


Combining this with (19), we get 


bl bl ^ r a / ■ t Cl • f Cl 

sup — = sup — < tan pi < ml — = mf —. 

n 0 ai n + ai n 0 h n + 6, 


Thus, the condition of Lemma |2.1| is satisfied, which implies that 

V • (Bi Vit) = <9 x (7i,o d x u ) + (pip <%«) + ^(71,2 ^u), 
where 71,0, 71,1, and 71,2 are nonnegative and given by 

(a — e) — (b + e/M) cot(/3i), for 6>0 


7i,o = 


for b < 0 


b + e/M 


1 R ■ R ’ for 6 > 0 
7!,! = ^ cos pi sin pi 

0, for b < 0 


7i,2 = 


(c—e) — (6 + e/M)tan/?i, for b> 0 
e, for b < 0. 


Similarly, for sufficiently large M, from © we have 


C2 


C2 


sup — = — M < sup — = sup 


c — e 


62 


62 


0+ b 2 
s ‘o 


b 2 


inf — = — — > inf — = inf-. 

q+ 02 M si- a2 a — e 


Combining this with (20) gives 


c 2 c 2 q ^ . , b 2 

sup — = sup — < tanp2 < mt — = mf —. 
a 0 b '2 q- b 2 n 0 a 2 n~ a 2 


( 20 ) 


( 21 ) 


Then from Lemma o we have 


V • (Da Vit) = ^(72,0 <9 x w) + dp 2 (72,1 dp 2 u) + d y ( 72,2 5 y tt), 


where 72,0, 72,1, and 72,2, all nonnegative, are given by 


72,0 = 


72,1 = 


J e + e/M cot/? 2 , 

I (a — e) — b cot fa 

' -e/M 

cos fa sin 02 ’ 
b 

, cos 02 sin fa ’ 


for b > 0 
, for b < 0 

for b > 0 

for b < 0 


( 22 ) 


6 















72,2 = 


e + e/M tan^2, 

(c — e) — b tan /?2, 


for b > 0 
for b < 0. 


By combining (18), ([21]), and ( |22| ), we get 


V • (BVu) = d x (~fo d x u) + 9 Pl u) + dp 2 u) 

where 70, 7+, 7J", and 72, all nonnegative, are given by 

+ dy{ 72 

d y u ), 

7 o = 7 i,o + 72,0 = | 

\a — (b + e/M ) cot + e/M cot @2, 

for 

b > 0 

for 

b < 0 

la — 6 cot /?2, 

72 = 71,2 + 72,2 = | 

\c — (b + e/M ) tan / 3 i + e/M tan/?2, 

for 

b > 0 

for 

b < 0 

7 i + = 7 i,i = 

1 c — 6tan/?2, 

i 7 e/ " fOT 620 

COS Pi Sin Pi 


0 , 


7i = 72,1 = < 


-e/M 


cos @2 sin /?2 
b 

cos /?2 sin fa 


for b < 0 
for b > 0 
for b < 0. 


Since the above result holds for any sufficiently small e, we take the limit as e —> 0 and obtain the 
conclusion of the lemma. □ 


In the next lemma, we show that (10) can be satisfied uniformly over Tt under a reasonable 
assumption on the diffusion matrix D. For a positive number R , we denote 


Br{xq, yfi) = {(x,y) (= £1 : \(x,y) — (x 0 ,y 0 )| < #}. 


Lemma 2.3. Assume that Q CM 2 is a bounded domain and D is continuous on Q = £1 U 
and symmetric and uniformly positive definite on fl. Then, there exists R > 0 such that for any 
(x 0 ,y 0 ) e Q, 


where 


K x iV) , « ^ f c(x,y) 
Bf(x 0 ,y 0 ) a ( X, y} ^a B+(x 0 ,y 0 ) b(x, y) 

c(x,y) a &(z 5 y) 

B~(xo,yo)b( x ’y) 3a B R (xo,yo) a \ x i V) 

b{x,y) , a ^ ■ t a(x,y) 

Bf(x 0 ,yo) C ( X ' ^ 3 a B+(x 0 ,yo) b\ x ,y) 

a( x ,y) , a b{x,y) 

sup u -f + mf -w 

B R (x 0 ,y 0 )^ X 'y> 3a B R (x 0 ,yo) c (x,y) 


(23) 

(24) 

(25) 

(26) 


a = inf (a(x, y)c(x, y) - b(x, yf) , 

a = max ( sup a(x, y)(\b(x, y)\ + 1), sup c(x, y)(\b(x, y)\ + 1) 
V q 


(27) 

(28) 
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Proof. The inequalities (25) and (26) can be obtained from (23) and (24) by the symmetry between 
a and c. To prove (23) and (24), we first notice that the uniform positive definiteness of D implies 
that a > 0 and 

a(x, y)c(x , y) - b(x, y ) 2 > a, V(.x, y) G Q. 

Dividing both sides by ab, we have 


a 


c(x,y) > b{x,y) _ 

b(x,y) ~ a(x,y) a(x, y)b{x, y)'' 
c(x,y) < b(x,y) a 

b{x,y ) “ a(x,y) a(x, y)b(x, y)'' 


for b(x, y) > 0 
for b(x, y) < 0. 


Denoting 


F(x,y) = 


c(x,y ) 
b(x, y) ’ 


G = 


from the definition of a (27) we have 


b(x,y ) 
a(x,y)' 


(29) 


OL 

F{x, y) > G(x, y) + 

a 

& 

F(x,y) < G(x, y) -, 

a 


for b(x, y) > 0 
for b(x, y) < 0. 


Define 


Of 

M = sup \ G(x, y)\ + 
n a 


(30) 


To prove (23), we introduce a cut-off function of F as 



'F(x, y), 

for 

b(x,y ) > 0 and F(x,y) < M 


F + (x,y) = < 

M, 

for 

b(x,y ) > 0 and F(x,y) > M 

(31) 


M, 

for 

b(x,y) < 0. 



This function is continuous on D and has the following properties, 

F + (x,y) < F(x,y), for 6(x, y) > 0 (32) 

77 / 

F + {x,y) > G(x,y) + -, for (x,y)€Q. (33) 

a 

Since D is bounded, F + and G are uniformly continuous on D. Thus, there exists a constant R\ > 0 
(which depends on F + , G, a, and a and therefore the diffusion matrix D but not on (xo,yo)) such 
that 

F + (x,y) > F + (x 0 ,y 0 ) ~ G(x, y) < G(x 0 , y 0 ) + V(x, y) G B Rl (x 0 , y 0 ) (34) 

6 a 6a 

which in turn implies that 


Q/ Q/ 

inf F + (x,y) > F + (x 0 ,y 0 ) - —, sup G(x, y) < G(x 0 , y 0 ) + 
B Rl (x 0 ,yo) oa B Rl (xg,yo) 


Combining these results with (33), we have 
sup 

B Rl (x 0 ,yo) 


2 a/ 

sup G(x,y) < G(xo,yo) + — < F + (x 0 ,yo) ~ — < inf F + (x,y)- 


6 a 


3a b r (x 0 ,yo) 


3 a 





























or 


a 


sup G(x,y)< inf F + (x,y)~ —. 

Br^^xqiUo) Br\\X Oj2/o) OQ? 


From this and (32), we have 


Ot 

sup G(x,y) < sup G(x,y) < inf F + (x,y)~ — 

Br.i(x o,yo) Bri ( x o >2/o) oCk 


B Rl {*o,yo) 


c\ a 

< inf F + {x,y )-< inf F(x,y) -, 

B+^x 0 , 2 / 0 ) 3 « S+^xo.yo) 3 « 


which gives (23) (with i? = R\). 
Similarly, we define 



F(x,y), 

for 

b(x, y) < 0 and F(x, y) > —M 


IS 

II 

-M, 

for 

b(x, y) < 0 and F(x, y) < —M 

(35) 



for 

d 

Al 

"S 

-0 



It is continuous on fl and satisfies 

F~(x,y) > F(x,y), forb(x,y)<0 


a 


F (x, y) < G(x, y) -, for (x, y) <5 fi. 

a 


(36) 

(37) 


From the uniform continuity of F and G on fi, there exists a constant R 2 > 0 (which depends 
only on the diffusion matrix D) such that 


Thus, 


_ _ a ol 

F (x, y) < F (xq, yo) + —, G(x,y) > G(x 0 ,y 0 ) - V(x,y) G B R 2 (x 0 ,yo)- 


o> 2 a a 

sup F~(x,y) < F~(xo,yo) + < G(x 0 ,y 0 ) - — < inf G(x,y)~ — 


(38) 


Br 2 (xo,2/o) 


3a 


3a B R „(xo,yo) 


3a 


or 


a 


sup F (x,y)< inf G(x,y)-—. 
B R2 (x 0 , yo ) Bfi 2 (*o,yo) 3a 


Then, the inequality (24) (with R = R 2 ) follows from this and (36). 


Using the symmetry between a and c we can show from (123) and (|24|) that (25) and (26) hold 


with some positive numbers R 3 and R 4 , respectively. Taking R = min(i?i, R 2 , R 3 , R 4 ), we have 
proven that (23) - (26) hold for the same R. □ 


Remark 2.2. The constants a and a defined in (27) and (28) and the functions F, F + , F , and 


G defined in (29), (30), (31), and (35) are determined completely by the diffusion matrix D. So 


is the constant R. If F + , F , and G are Lipschitz continuous with constants L F +, L F ~, and Lq , 
respectively, then we can choose R as 


R = 


a 


3a max(L F +, L F ~. Lq) 


(39) 

□ 
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Theorem 2.1. Assume that 12 C M 2 is a bounded domain and the diffusion matrix D is con¬ 
tinuous on 12 and symmetric and uniformly positive definite on 12. Then, there exists a constant 
R > 0 (which depends only on DJ such that for any (xo>2/o) £ the intervals 


sup 


b(x,y) 


inf 


c{x,y) 




{xo , yo) B+(x 0 ,y 0 )b{x,y) 


sup 

\B~(x 0 ,yo) 


c{x,y) . nf b(x,y) 
K x iVY B-(x 0 , yo )a(x,y) 


are not empty when B R (xo,yo) and B R (xo,yo) are not empty, respectively. Moreover, for any 
constants fii and @2 with tan /3\ and tan @2 in the intervals, the nonnegative directional splitting 
0 holds with 12q = Bji(xo,yo) and nonnegative coefficients given in (Tify - |l5|). 


Proof. The theorem is obtained by combining Lemma |2.2| and Lemma 2.3 


□ 


Remark 2.3. The result in the above theorem is different from Weickert’s result (31] in the 
sense that the coefficients in the directional splitting ©> are nonnegative in a neighborhood of 
(* 0 ) 2 /o) whereas those of Q are nonnegative only at (xo,yo)- This is crucial to the construction 
of monotone FD schemes since any FD scheme for the divergence form uses the values of the 
coefficients on a neighborhood or a stencil of a mesh point for the discretization at the point. □ 


Remark 2.4. As mentioned in Remark 2.2 the constant R depends on the diffusion matrix D 
and is independent of the location. This implies that a uniform mesh can be used to construct a 
monotone FD scheme as long as the mesh spacing is sufficiently small and the maximal size of the 
FD stencils (used for all mesh points) is bounded (see ^3l). □ 


3 Construction of a monotone FD discretization 


In this section we develop a monotone FD discretization for ([Tj) based on the nonnegative directional 
splitting pTj ). We consider a uniform mesh {xj,yf) = ( jh , kh ), j, k = 0,..., N, where N is a positive 
integer and h = 1/N. 

Let ( Xj,yk ) be an interior mesh point. Denote by ca FD stencil of size (277^ + 1) x (2m J j c + l) 
centered at ( Xj,yk ). We assume that 


m jjk < 


3a 

a 


+ 1 ) 


(40) 


where a and a are defined in (27) and (28) and [x] denotes the nearest integer smaller than or 
equal to x. We will show that the right-hand side of ( |40[ ) is an upper bound of the stencil size 
to guarantee a monotone FD discretization for ([!]). We also assume that the mesh spacing h is 
sufficiently small such that the stencil ojjj, C Bfi(xj,yk), with R being the radius determined in 
Lemma |2.3| and Theorem |2.1[ Mathematically, this is equivalent to 


\f 7 lh max rrto fc < R. 
j,k 


Using (39) and (40), we obtain a sufficient condition for (41) as 


V2h 


3a 


a 


+ 1 < 


a 


3a max(L^+, Lp -, Lq) 


(41) 


( 42 ) 


which can be achieved by choosing h sufficiently small. 
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We first consider the case where Uj k C 0. Under the above assumptions, Theorem 2.1 implies 
that the nonnegative directional splitting 0 , with the coefficients being nonnegative on Uj ik , holds 
for any constants j3\ and fa satisfying 


_ b(x,y) _ c(x, y) 

Aj k = sup — -- < tan pi < B j k = mf —- 

.,+ a(x,y) u + b[x, y) 


i,k 


_ c(x, y) _ b(x, y) 

C— sup —-- < tan p 2 < Vj k = mi 


j,k 


b(x,y) 


■~ k a(x,y) 


(43) 

(44) 


It is not difficult to see that 
ci(x,y) 


A 7i = inf 


o-i b(x, y ) 

B -k =sup 


C~l = inf 


+ b[x, y) ’ j ’ k c(x, y) 


j,k 




From Lemma 2.3, we have 


Bj,k 

Dj,k 


' / 

c{x 

,y) 


a 

k ^ 



3 a 


a 

,fc > 



b(x,y ) ! a(x,y) 

D- k = sup 


j,k uj k c(x,y)’ j ’ k J b{x,y) 

7 ,k 


A~ 1 - B _1 > — 

A j,k *j,k A 3a , 


3a 


C~} - D~l > —, 

3, k 3,k ~ 3a ’ 


j,k 


(45) 

(46) 


where a and a are defined in (27) and (28) 


The central FD discretization of the first and fourth terms of © is straightforward. It can be 
done based on points (xj- 1 } y k ), (xj,y k ), (x j+1 ,y k ) and (xj, y k -i), ( Xj,y k ), (xj,y k +i), respectively. 
(These points are on a 3 x 3 stencil centered at [xj. y k ) ) On the other hand, the FD discretization 
of the second and third terms are more complicated. Take the second term as an example. The 
main challenge of the discretization is to find fi\ with tan/3j 6 ( Aj k ,Bj k ) such that the straight 
line having the slope tan/?i and passing the point ( Xj,y k ) intersects at least two more mesh points 
on ujj jk other than ( Xj,y k ). Once such a /3 1 has been found, the three mesh points are then used 
to discretize the term (which is the second order directional derivative along tan/?i) using central 
finite differences. 



Figure 1: Sketch of a FD stencil and principal directions. 

To find /?i, we follow [31] to define 4 mj jk principal directions associated with the outer mesh 
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nodes of ujj. k (cf. Fig. [I]) 


A = 


arctan 

) , 

for 


J 


arctan 

' m i,k \ 

for 

1 2 J ’ 

arctan 

' m .i,k \ 

for 

y-2 m jtk -i) > 


(47) 


- 1. 


When mj } k is sufficiently large, we can choose Al to be one of A; i = 1, 2rrij_ k that satisfies (43). 
Indeed, there are three possibilities Aj k and Bj k are located: Aj k < 1 < Bj k , Aj k < B jk < 1, or 
1 < Aj,k < Bj t k- For the case Aj k < 1 < Bj tk , we can choose /?i = f3 mj k = 7t/ 4 and rrij^ = 1- It is 
easy to see that (43) is satisfied and the three mesh points on the diagonal line running from the 


southwest corner to the northeast corner can be used to discretize the directional derivative. 

For the case Aj k < Bj ik < 1, we need to consider the principal directions A> i = 1,..., rrij^ — 1. 
Thus, we require that at least there is an i (1 < i < rrij tk — 1) such that 


A j-k < 


k^j.k 


< B 


j ,k- 


(In this case, Al is chosen as A-) Equivalently, nij yk should be sufficiently large such that (mj jk Aj k , rrij^Bj ^) 
contains at least an integer. A sufficient condition for this is that the length of the interval is bigger 
than one, i.e., 

m i,k > ~B - A • ( 48 ) 

&j,k — A j,k 

For the case 1 < Aj >k < Bj jk , we need to consider the principal directions i = rrij tk + 

1,..., 2mj^ — 1. We require 

A ^ m j’ k ^ p 

A i’ k < o.- „ < B i’ k 


2 m jt k - i 


or 


m j, kB j t k < 2 m j,k - i < m j,kAjJ, 


-1 


for some rrij^ < i < 2mj )k . This is mathematically equivalent to the requirement that (m Ak B rj k , m J:k A- ' 
contains at least an integer. A sufficient condition for this is 


mj,k > A - 1 - B - 1 ' 

■^jjk n j,k 

The above analysis shows that we should choose mj )k such that 


Wlj,k — 1 ) 

(mj t kAj^,rrij t kBj t k) contains at least an integer, 
^ (rnj kBj^. mj kAj^) contains at least an integer, 

A sufficient condition for this is 


for Aj, k < 1 < Bj :k 
foi A jj ; <t Bj _ k ' ■ 1 
for 1 < A ]k < Bj k- 


(49) 


' m j,k > 1, 

for A jtk < 1 < B j k 

m 3A > B ]ik -A jik > 

for A jtk < B jt k < 1 

m j, k A A -i B -i, 

3,k D j,k 

for 1 < A jt k < B j k 


( 50 ) 
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From (45), one can see that the above condition can be satisfied when rrij^ is taken as rrij tk = 
[lr] + 1- Since this is only a sufficient condition, so we have (40). 

Similarly, the condition that we can choose @2 as one of the principal directions /%, i = 
— 2mj t k, — 1 satisfying (44) is 


Tflj.k — 1 ) 

< {rrij t kCj t k-, m j,kDj ) k) contains at least an integer, 
X m j,kDjk-,mj^Cjk) contains at least an integer, 

A sufficient condition for this is 


for C jtk < -1 < D jtk 
for 1 f'y .k ^ D 

for Cj :k < D j k < -1. 


m^k > 1, 

for Cj^k ^ 1 -^j,k 

m j,k > D jik -C jik ’ 

for 1 ^ ^ Dj^k 

> ~n~T D -i ) 
K -'j,k j,k 

T—1 

1 

VI 

-g 

Cl 

V 

-se 

a 


(51) 


(52) 


From (46), one can see that the above condition can be satisfied when rrij )k is taken as the right-hand 
side ofl40b. 



Figure 2: FD stencil near the domain boundary 


We now consider the situation when the stencil ojj k is not in Q. This can happen when ( Xj,y k ) 
is close to the domain boundary and 'mj.k > 1. In this case, we can determine rrij}.. fii, and /?2 as 
for the case where u)j k is in Q. But, one of the three mesh points found along direction tan j3\ or 
tan /3'2 lays outside of the domain; see Fig. [2] A treatment for this difficulty is to replace the mesh 
point by the intersection of the straight line (with slope tan f3\ or tan fif) and the domain boundary 
and then discretize the corresponding directional derivative using central finite differences based 
on this boundary point and other two mesh points. 

It is worth pointing out that each of the terms on the right-hand side of © has been discretized 
using central three-point finite differences along the corresponding direction. Since the coefficients of 
the differential operator are nonnegative, it is not difficult to see that the coefficient matrix of the FD 
equations is a Z-matrix (with off-diagonal entries being nonpositive) and has nonnegative diagonal 
entries. Moreover, it is easy to verify that the matrix is diagonally dominant and irreducible. Thus, 
the coefficient matrix is an M- matrix and the FD scheme is monotone. 

For easy reference, we summarize the above analysis in the following theorem. 

Theorem 3.1. Assume that the diffusion matrix D is continuous on D and symmetric and 
uniformly positive definite on D and a uniform mesh with spacing h is given for D. Ifhis sufficiently 
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small, a monotone FD scheme exists for problem & and can be constructed based on central 
(directional) finite differences. 


Remark 3.1. Generally speaking, the upper bound for m 7 -fc defined in (40) is very conservative. 


In practice, we can use stencils of variable size for different mesh points. The stencil size mj,h can 


be found using (|49|) and (|5T|) . The constant R can be calculated using (39). The mesh spacing can 

(53) 


be verified using 


\ff.h max m, j, < R. 
j,k 


□ 


To conclude this section, we would like to see what condition is if only the three-by-three 
(jrij k = 1) stencil is used for all interior mesh points. In this case, we need to use the mesh points 
on the diagonal lines of ojj % k to discretize the second and third terms on the right-hand side of m. 
This means that Ajk < 1 < Bjk (tan/3i = 1) and Cjj. < —1 < Dj^ (tan/32 = —1). From (43) and 
(44), this implies 

a(x,y) > \b(x,y)\, c(x,y) > \b(x,y)\, V(i,//)eO 


or, in words, 
theorem. 


is strictly diagonally dominant for any point in II. 


(54) 

Thus, we have the following 


Theorem 3.2. Assume that the diffusion matrix D is continuous on D and symmetric and 
uniformly positive definite on It and a uniform mesh with spacing h is given thereon. If {51)) is 
satisfied and h is sufficiently small, then a monotone FD scheme can be constructed for & based 
on a three-by-three stencil and central finite differences. 


It is pointed out that the condition (54) has been obtained by Weickert m and Mrazek and 
Navara 


4 Numerical results 

In this section we present some numerical results for four examples to verify the theoretical findings 
in the previous sections. We first note that it is useful to have an estimate on the values of 


parameters R (the uniform radius defined in Lemma 2.3) and rrijj, (the size of FD stencil) in the 
actual computation. The former can be estimated using (39) where a, a, L F +, Lp~, and Lq can 


be approximated numerically. For the latter, (40) is generally too conservative to use in practical 
computation. Instead, we use (49) and (51) to estimate with A 3 y., BjCj^, Dj^ being 

computed on Bp(xj,yk). Notice that we use Bp(xj,yk) instead of ujj t k as defined in (|43|) and 


This is because the former is larger than the latter and the former is independent of rrijk which is 
to be sought. Once mj,k has been determined, we can find two principal directions (cf. ©) to 
satisfy (43) and (44) and finally construct the FD discretization of 0- 

Example 4.1. This example is in the form of ([Tj) with 


9 4sin(27 xxy) 

4sin(27 rxy) 3 

We seek to verify that the developed monotone scheme satisfy DMP and choose 


(55) 


f(x,y) = 0, g(x, y) = cos(nxy) + y. 


14 






























This example does not satisfy (54). So we do not expect that a 3 X 3 stencil will work for this 
example. A direct numerical calculation shows that 


a = 11 , 


a = 45. 


Thus the upper bound in (40) is [2 a/a] + 1 = 13, which requires a stencil of size 27 x 27 for the 
monotone FD discretization. On the other hand, a numerical calculation based on ( |49[ ) and © 
shows that a stencil of size 5 x 5 is sufficient to guarantee a monotone FD discretization. 

The numerical results obtained with a stencil of size 5x5 for this example are shown in Table [l] 
They confirm that the computed solution satisfies the maximum principle and stays between the 
minimum and maximum values of the solution on the domain boundary. 


J = I< 

Boundary Min 

Interior Min 

Boundary Max 

Interior Max 

21 

-5.105652e-2 

1.040961e-2 

2.000000 

1.912261 

51 

-5.105652e-2 

-2.444841e-2 

2.000000 

1.972163 

101 

-5.105652e-2 

-3.753179e-2 

2.000000 

1.987642 

501 

-5.109830e-2 

-4.837255e-2 

2.000000 

1.997862 

1001 

-5.110190e-2 

-4.973654e-2 

2.000000 

1.998961 


Table 1: Example 4.1 Extrema of the numerical solutions in the interior and on the boundary of 
the physical domain. 


Example 4.2. In this example, we seek to verify the accuracy of the monotone scheme. We use 

in the form of ([Tj) with (55), but choose / 


the same differential operator as in Example 4.1 


i.e., 


and g such that the exact solution of the BVP is given by u(x, y ) = sin(27rx) sin(37n/). 

The numerical results obtained for this example are shown in Figs. [3] and [4j The contours of 
the exact and computed solutions are shown in Fig. [3j The sign pattern in Fig. |4a| shows that the 
coefficient matrix has nonpositive off-diagonal entries and nonnegative diagonal entries, as stated 
in Theorem 3.1 Since the scheme is monotone, from the result of Barles and Souganidis [jj we 


may expect the computed solution to converge in second order, which indeed can be seen from the 


convergence history shown in Fig. 4b 


Example 4.3. For this example, the diffusion matrix is given by 

1.1 sm.{2itxy) 
sin(27rx?/) 1.1 


(56) 


The functions / and g are chosen such that the exact solution of the BVP reads as u(x , y) = 
sin(27rx) sin(37ry). This example satisfies (54) and by Theorem |3.2[ we can use a 3 x 3 stencil for 
constructing a monotone scheme. 

The numerical results for this example are shown in Figs. [5] and |6j One can see that the 
coefficient matrix has nonpositive off-diagonal entries and nonnegative diagonal entries and the 
computed solution converges to the exact solution at a rate of second order. 

Example 4.4. In this final example, the diffusion matrix is given by 


cos 6 — sin 6 


k 0 


cos 6 sin 6 

sin 6 cos 6 


0 1 


— sin 6 cos 9 


(57) 


where 6 = 7rsin(x) cos (y) and /c is a parameter. The bigger k is, the more anisotropic the diffusion 
matrix is. The more anisotropic D is, the larger stencil is required for the construction of a 
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(a) Exact solution 


Figure 3: Example 


4.2 


(b) Numerical solution with J = 81 
Contours of the exact and computed solutions. 




(a) Sign pattern for J = 81 (b) Convergence history 


Figure 4: Example 4.2 Sign pattern of the coefficient matrix and the convergence history of the 


maximum norm of the error of the FD solution. 
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(a) Exact solution 


Figure 5: Example 


4.3 


(b) Numerical solution with J = 51 
Contours of the exact and computed solutions. 



nz = 2601 

(a) Sign pattern for J = 51 



(b) Convergence history 


Figure 6: Example 4.3 Sign pattern of the coefficient matrix and the convergence history for the 
maximum norm of the error of the FD solution. 
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Figure 7: 
solution. 


Example 


4.4 


The convergence history for the maximum norm of the error of the FD 


monotone FD scheme. We consider two cases with k = 10 and k = 100. For both of these cases, 
we again choose / and g in ([Tj) such that the exact solution is u(x,y) = sin(27ra:) sin(3ny). 

The numerical calculation suggests that a 7 x 7 stencil be needed to construct a monotone 
scheme for k = 10 and a stencil of 53 x 53 for the case k = 100. The second order convergence can 
be seen in Fig. [7] for both cases although the error is almost two magnitude larger for the much 
more anisotropic case with k = 100. 


5 Conclusions 


In the previous sections we have studied the construction of monotone finite difference schemes 
for the divergence form (TO) using nonnegative directional splittings. The main results are stated 

Theorem 2.1 states that the diffusion operator in ([Tj) has a non¬ 


in Theorems 2.1, 3.1 and 3.2 


negative directional splitting at any arbitrary interior point (cf. ©) when the diffusion matrix 
D is continuous on Q and symmetric and uniformly positive definite in Q. This splitting holds 
in a neighborhood of the point in the sense that the corresponding coefficients are nonnegative in 
the neighborhood. Moreover, the size of the neighborhood depends only on D. This result is an 
important improvement over that in m where the coefficients of the splitting are nonnegative only 
at a point. Another improvement is that the function b(x,y) in Q is allowed to change sign on 
D. These improvements are crucial to the construction of monotone finite difference schemes for 
([!]) since any finite difference discretization for the divergence form (JTj) requires information of the 
coefficients in a neighborhood of a mesh point. 

The construction of monotone finite difference schemes using the nonnegative directional split¬ 
ting 0 has been discussed in l[3j A key to the construction is to define principal directions for the 
underlying finite difference stencil m ■ We have shown that some of these principal directions can 
be used as the splitting directions when the mesh spacing is sufficiently small and the size of finite 
difference stencil is sufficiently large. Note that the stencil size has a uniform bound determined 
by D; see (40). This bound can become large when D is strongly anisotropic. In that case, not 
only a large stencil but a very fine mesh have to be used. In such a situation, the constructed 
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finite difference scheme may be impractical and a more sophisticated way to construct monotone 
finite difference schemes may be needed, which is a research topic worth more investigations in 
the future. On the other hand, Theorem 3.2 shows that a 3 x 3 stencil can be used to construct 
monotone finite difference schemes when D is strictly diagonally dominant. Numerical results in ^4] 
are consistent with the theoretical findings. 
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